Анализ результатов моделирование самосборки липидного бислоя

Автор: Федоров А.Н.

In [1]:
import sys
import pymol

_stdouterr = sys.stdout, sys.stderr
pymol.finish_launching(['pymol', '-q'])
sys.stdout, sys.stderr = _stdouterr

from pymol import cmd, movie
In [2]:
import time
from IPython.display import Image, HTML


def render(width=640, height=480, filename="/tmp/tmp.png", sleep=1):
    cmd.ray(width, height)
    cmd.png(filename)
    time.sleep(sleep)
    return filename

No module named 'Pmw'
Unable to initialize plugin 'apbs_tools' (pmg_tk.startup.apbs_tools).

Подготовка исходных данных

Загрузим все предоставленные в инструкции файлы в текущую папку:

In [1]:
files = ['dppc.itp', 'lipid.itp', 'dppc.gro', 'b.top', 'em.mdp', 'pr.mdp', 'md.mdp']
for f in files:
    !wget -q http://kodomo.cmm.msu.ru/~golovin/bilayer/{f}

На основе одного липида созадим ячейку с 64 липидами:

In [ ]:
!genconf -f dppc.gro -o b_64.gro -nbox 4 4 4

Посмотрим на наш липид и финальную ячейку в pymol:

In [ ]:
%%bash
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
In [8]:
cmd.reinitialize()
cmd.load('dppc.pdb')
Image(render())
Out[8]:
In [10]:
cmd.reinitialize()
cmd.load('b_64.pdb')
Image(render())
Out[10]:

Все хорошо, результаты ожидаемы.

В текстовом редакторе в файле b.top установим правильное количество липидов в системе:
// 32 -> 64

Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды:

In [ ]:
!editconf -f b_64.gro -o b_ec -d 0.5

Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул:

In [ ]:
%%bash
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v

Начальное значение максимальной силы: Fmax=437970.

Конечное значение максимальной силы: Fmax=616.887.

It works, максимальная сила уменьшилась на порядки.

Добавим в ячейку молекулы воды типа spc:

In [ ]:
!genbox -cp b_em -p b -cs spc216 -o b_s

Проведём "утряску" воды:

In [ ]:
%%bash
grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
In [ ]:
%%bash
grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v
In [ ]:
%%bash
grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v

Посмотрим на b_pr и b_s в pymol:

In [ ]:
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
In [13]:
cmd.reinitialize()
cmd.load('b_s.pdb')
Image(render())
Out[13]:
In [14]:
cmd.reinitialize()
cmd.load('b_pr.pdb')
Image(render())
Out[14]:

Нельзя сказать, что здесь заметны какие-то сильные отличия. Все +- похоже, лишь появилось чуть больше молекул воды между липидами.

Запуск моделирования

Для начала протестируем, что все работает локально:

In [ ]:
!grompp -f md -c b_pr -p b -o b_md -maxwarn 1
In [ ]:
!mdrun -testverlet -deffnm b_md -v -maxh 0.002

Вроде бы все хорошо. Делаем заявку через google-doc.
Path: /home/anfedorov/hw7
tpr file: b_md.tpr

Иии... Ждем.

Визуализация динамики системы

Подготовим анимацию системы:

In [ ]:
!trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
In [23]:
!mkdir -p ./frames/
!rm ./frames/*

cmd.reset()
cmd.reinitialize()
cmd.load("b_pbc_1.pdb")
cmd.mpng('./frames/movie')
In [ ]:
!ffmpeg  -framerate 24 -y -i './frames/movie%04d.png' -c:v libx264 -pix_fmt yuv420p -vf "crop=trunc(iw/2)*2:trunc(ih/2)*2" movie.mp4 > /dev/null
!rm ./frames/*
In [26]:
import base64

with open("movie.mp4", 'rb') as file:
    video = base64.b64encode(file.read()).decode()

video = '<video controls alt="test" src="data:video/mp4;base64,{}">'.format(str(video))
HTML(data=video)
Out[26]:

Исходная сетка разрушается сразу, бислой начинает формироваться с ~10 модели (t = 4500) и заканчивает к ~25 модели (t = 12000).

Площадь одного липида

In [ ]:
!g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg

Нормалью к бислою является ось Х. Подсчитаем изменение числа липидов на единицу площади бислоя:

In [5]:
import matplotlib.pyplot as plt
import pandas as pd


df = pd.read_csv('box_1.xvg', header = None, skiprows = 24, sep='\t')
T, X, Y, Z = df[1], df[2], df[3], df[4]

nlipids = 64
plt.plot(T, Y * Z / nlipids)
plt.gca().set(xlabel="Time", ylabel="$\\frac{S_{YZ}}{N_{lipids}}$")
plt.show()

Ожидаемо, число липидов на единицу площади ячейки падает - происходит компактизация системы, липиды организуются в искомый бислой.

Гидрофобная и гидрофильная поверхность

In [ ]:
!g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100

Построим зависимость доступной растворителю гидофильной/гидрофобной поверхностей от времени:

In [65]:
df = pd.read_csv('sas_b.xvg', header = None, skiprows = 22, sep='\s+')
T, phobic, philic = df[0], df[1], df[2]

plt.plot(T, phobic, label='Hydrophobic')
plt.plot(T, philic, label='Hydrophilic')
plt.gca().set(xlabel="Time", ylabel="Solvent Accessible Surface (nm\\S2\\N)")
plt.legend()
plt.show()

Теоретически, липидный бислой является более энергетически выгодной структурой, чем просто случайное расположение молекул. Это связано с минимизацией энергетически невыгодных контактов молекул воды и гидрофобных хвостов липидов.

Наши теоретические идеи подтверждаются как качественно(визуализация), так и количественно(график). Видно, что при сборке бислоя резко падает площадь доступной растворителю гидрофобной поверхности. Гидрофобные хвосты липидов группируются вместе, минимизируя число контактов с водой.

Оценка фазового состояния

Скачаем индекс и запустим gromacs:

In [ ]:
!wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
!g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
!g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X

Строим зависимость:

In [66]:
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')

plt.plot(start[0], start[3], marker='o', label='Start')
plt.plot(end[0], end[3],marker='o', label='End')
plt.gca().set(xlabel="atom", ylabel="oder measure")
plt.legend()
plt.show()

В финальном бислое атомы значительно более упорядочены чем в начале траектории. Это ожидаемо, т.к. бислой - это упорядоченная структура с характерным расположением липидов.

Выводы

Получилась довольно интересная работа, наглядно показывающая сборку бислоя, фундамента клеточных стенок. Интересно, можно ли смоделировать не только небольшой участок бислоя, но, например, комплекс бислой + мембранный белок и сравнить с готовой pdb моделью (белок, для простоты моделирования, можно брать в уже нативной конформации).

Теоретически, это все еще будет энергетически выгодно, т.к. мембранных белки имеют гидрофобные участки, которые и интегрируются частично/полностью в липидный бислой. А результаты моделирования наглядно бы показали, насколько мы можем доверять МД в подобных задачах.

In [68]:
!jupyter nbconvert hw7.ipynb --to html

[NbConvertApp] Converting notebook hw7.ipynb to html
[NbConvertApp] Writing 5108588 bytes to hw7.html
In [ ]: